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The growth of a well-formed epithelial structure is 
governed by mechanical constraints, cellular apico- 
basal polarity, and spatially controlled cell division. 
Here we compared the predictions of a mathematical 
model of epithelial growth with the morphological analy- 
sis of 3D epithelial structures. In both in vitro cyst models 
and in developing epithelial structures in vivo, epithelial 
growth could take place close to or far from mechani- 
cal equilibrium, and was determined by the hierarchy 
of time-scales of cell division, cell-cell rearrangements, 
and lumen dynamics. Equilibrium properties could be 



inferred by the analysis of cell-cell contact topologies, 
and the nonequilibrium phenotype was altered by inhibit- 
ing ROCK activity. The occurrence of an aberrant multilu- 
men phenotype was linked to fast nonequilibrium growth, 
even when geometric control of cell division was correctly 
enforced. We predicted and verified experimentally that 
slowing down cell division partially rescued a multilumen 
phenotype induced by altered polarity. These results im- 
prove our understanding of the development of epithelial 
organs and, ultimately, of carcinogenesis. 



Introduction 

Epithelial morphogenesis is a complex process involving cell 
divisions, cell-cell and cell-ECM adhesion, cell migration, cell 
shape changes, and apoptosis, and represents a fundamental step 
in organogenesis. Indeed, these features are fundamental for the 
correct functioning of the tissue in terms of proliferation, sur- 
vival, and differentiation. Aberrant epithelial architecture is most 
frequently found in the pathogenesis of epithelial tumors, and 
architectural patterns have been used for decades by pathologists 
to diagnose and classify carcinomas. The study of morphoge- 
netic processes leading to the formation of epithelial tissues can 
thus be used to gain a better understanding of the development of 
epithelial organs and of carcinogenesis. 
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In vitro biological models have been successfully used to 
reproduce some of the key events involved in epithelial morpho- 
genesis, and represent a fundamental tool to dissect the molecu- 
lar cascade of events leading to the formation of tissues (O'Brien 
et al., 2002; Debnath and Brugge, 2005). 

Cystogenesis is one of the best studied examples of epi- 
thelial morphogenesis in vitro (McAteer et al., 1986; O'Brien 
et al., 2002) and is considered to be a prototype for the develop- 
ment of several spherical structures encountered in vivo, such as 
acini, follicles, ampullae, and alveoli. Cysts are spherical mono- 
layers of epithelial cells enclosing a central lumen (McAteer 
et al., 1986). Cells within cysts are connected by specialized 
junctions and cell-cell adhesion structures lying in the basolateral 
sides, whereas a strong apicobasal polarization characterizes the 
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external surface, contacting the ECM, and the apical surface, 
facing the lumen. The correct architecture and the formation and 
maintenance of the lumen are crucial for normal cyst morphol- 
ogy and are altered in several common human diseases such 
as polycystic kidney disease (Boletta and Germino, 2003), hy- 
pertension (Iruela-Arispe and Davis, 2009), and many epithelial 
cancers, such as prostate carcinomas or preinvasive epithelial 
lesions (Debnath and Brugge, 2005). 

Despite the specificity inherent to diverse types of tissues, 
recent findings support the idea that the formation of several 
spheroidal epithelial structures could be generated by common 
mechanisms, and that shared features underlie the appearance 
of aberrant phenotypes (Datta et al., 201 1). 

The first general key aspect involved in the process of cyst 
growth is the mechanics of cell contacts. Epithelial cells are 
physically connected to the ECM via integrin receptors (O'Brien 
et al., 2001), and neighboring cells are tightly linked by cell-cell 
junctions via adhesion receptors, such as cadherins and nectins 
(Harris and Tepass, 2010). Cell shape variations are caused by 
local deformations of the cortical actomyosin network. The cu- 
mulative effect of differential cell-matrix and cell-cell adhe- 
sion processes and of cortical elasticity can be described in 
terms of interfacial tensions, which have been shown to be 
the driving force behind tissue formation in several biological 
models (Kafer et al., 2007; Lecuit and Lenne, 2007; Manning 
et al., 2010). 

A second aspect involves apico-basal polarization and the 
de novo generation of a luminal space. Luminogenesis proceeds 
through a coordinated series of molecular events starting with 
the exocytosis of apical membrane proteins (such as Crumbs3a 
[Crb3], podocalyxin [PCX], and Mucin 1 [Mucl])tothecell sur- 
face, leading to the formation of the nascent lumen in a region 
termed the apical membrane initiation site (AMIS; Schliiter 
et al., 2009; Bryant et al., 2010). Similar structures have been 
observed during vascular lumen formation in developing mouse 
aorta (Strilic et al., 2009) and during neural rod formation in 
zebrafish (Tawk et al., 2007). After the formation of the AMIS, 
an asymmetric distribution of the phosphoinositides PIP2 and 
PIP3 is established (Shewan et al. , 20 1 1 ). In particular, the apical 
region is enriched in PIP2 and PTEN, whereas PIP3 is localized 
exclusively to the basolateral membrane. The AMIS matures to 
form a preapical patch (PAP), and eventually a lumen expands 
(Martm-Belmonte et al., 2007; Ferrari et al., 2008; Bryant et al., 
2010; Datta etal., 2011). 

A third aspect is the spatial control of cell division. The 
apico-basal polarization of specialized molecules such as PIP2, 
PTEN (Martm-Belmonte et al., 2007), Cdc42 (Jaffe et al., 2008), 
the Cdc42-specific exchange factors Tuba (Qin et al., 2010) and 
Intersectin-2 (Rodriguez-Fraticelli et al., 2010), Par3 (Hao et al., 
2010), aPKC (Qin et al., 2010), and LGN (Zheng et al., 2010) 
restricts the formation of the mitotic spindle to the monolayer 
plane, ensuring the stability of the polarized architecture during 
subsequent cell divisions. 

Remarkable progress in understanding the molecular 
cascade giving rise to the formation of a cyst has been made 
thanks to 3D cultures of MDCK cells and other similarly well- 
polarized epithelial cell lines, which provide an excellent model 



for epithelial morphogenesis in vitro and give valuable insights 
into the basic processes driving the formation of the correct 
tubular architecture and its disruption in cancer (Debnath and 
Brugge, 2005; Datta et al., 2011; Lo et al., 2012). Nonetheless, 
the precise connection between molecular and mechanical cues 
and the resulting dynamics is still largely unknown. In particu- 
lar, it has not been demonstrated whether the simple aforemen- 
tioned mechanistic principles are sufficient to reproduce cyst 
architectures. Furthermore, although it is clear that misorienta- 
tion of the spindle is an ingredient leading to the aberrant mul- 
tilumen phenotype (Zheng et al., 2010; Engelberg et al., 201 1), 
it is not clear how and to what extent cyst architecture and ori- 
entation of cell division are related. 

To shed light on these issues, we developed a mathemati- 
cal model based on the phenomenology of cell-cell and cell- 
matrix interactions that reproduces the experimentally observed 
phenotypes in both the normal and the aberrant cases. Our 
model indicates that cystogenesis occurs through a sequence of 
states that are far away from mechanical equilibrium, and pro- 
vides evidence that the appearance of the aberrant multilumen 
phenotype is determined by misorientation of the spindle only 
in conjunction with strongly out-of-equilibrium dynamics. We 
verify our predictions by performing computer-assisted, full 3D 
reconstruction at different stages of growth of MDCK cysts 
and of 3D epithelial structures derived from mouse embryos. 
We demonstrate that both cysts in vitro and developing epithe- 
lial structures in vivo share the same nonequilibrium phenotype. 
Our model predicts that equilibrium can be favored by slowing 
down cell division, by altering cell-cell and cell-ECM inter- 
actions and cortical contractility, or by increasing cell motility. 
Moreover, we demonstrate experimentally that slowing down 
cell division partially rescues a multilumen phenotype, thus 
confirming our predictions. We provide experimental evidence 
that inhibition of ROCK activity drives the appearance of equi- 
librium configurations in MDCK cysts. Finally, our model sug- 
gests that an increase in cell duplication rates, a typical hallmark 
of hyperproliferative diseases like cancer, could alone be re- 
sponsible for the appearance of aberrant phenotypes, indepen- 
dent of the control of spindle positioning. 

Results 

Theoretical modelling of cyst: 
morphogenesis predicts a crucial role for 
cell duplication and mechanical relaxation 

The geometric structure of epithelia can be described in terms of 
the mechanical forces that result from cell-cell and cell-ECM in- 
teractions. The action of the actomyosin cortex can be described 
by an overall surface tension that acts so as to minimize the total 
cell surface (Lecuit and Lenne, 2007; Manning et al., 2010). 
Conversely, basal and lateral adhesion receptors favor the in- 
crease of contact surfaces (Kafer et al., 2007; Lecuit and Lenne, 
2007; Manning et al., 2010). These ideas have been successfully 
used to quantitatively describe the geometrical configuration of 
ommatidia (Hilgenfeldt et al., 2008) and the growth of wings 
(Classen et al., 2005; Hufnagel et al., 2007; Lecuit and Lenne, 
2007; Farhadifar et al., 2007) in Drosophila melanogaster and 
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the spreading of murine sarcoma cell aggregates in vitro (Douezan 
etal.,2011). 

We develop a simple theoretical model of a growing cyst 
by extending the reasoning presented in Manning et al. (2010) to 
a 3D context. Such a model, under the assumption of equilibrium 
dynamics, predicts that the observed configurations of a growing 
cyst are those that minimize mechanical energy (see Materials 
and methods and Eqs. 1 and 2) containing contributions from 
both elastic and adhesive processes. Indeed, minimizing this en- 
ergy under the assumption of spherical symmetry, it is possible to 
find analytically the equilibrium configuration for a cyst-like cell 
aggregate of any number of cells (see Materials and methods). 
This equilibrium configuration has a central lumen and a soccer 
ball-like structure of the kind described in Cox and Flikkema 
(2010). Under more general conditions, an analytical solution is 
impossible, and it is necessary to resort to computer simulations. 
To this aim, we implemented a dynamical Potts model (Potts, 
1952; Wu, 1982; Graner and Glazier, 1992; Glazier and Graner, 
1993) with energy given by Eqs. 1 and 2. Within this model, ran- 
dom motility of cells plays a role analogous to temperature, and 
the associated thermal energy is always much smaller than the 
typical values of mechanical energy. (A more detailed descrip- 
tion of the computational model is given in the Materials and 
methods section.) 

The dynamics governed by Eqs. 1 and 2 depend crucially 
on the mechanical relaxation timescale T, i.e., on the time 
needed to reach global equilibrium, and on the characteristic 
duplication time t. Whenever the ratio x = Tfr of the relaxation 
time to the division time is much smaller than unity (x ^ 1), 
cystogenesis proceeds through a sequence of equilibrium states. 
Conversely, when x is much larger than unity (x ^» 1), the next 
division round takes place before significant mechanical adjust- 
ments occur. In the latter case, the dynamics are constantly kept 
out of equilibrium by subsequent cell divisions. A deeper under- 
standing of the problem comes from the observation that con- 
strained systems often display a spectrum of relaxation times 
associated with a complex energy landscape with multiple local 
minima (Mezard et al., 1987). For the present purposes it will be 
sufficient to distinguish between the fast timescales that govern 
mechanical rearrangements without topology changes, such as 
the relaxation of the basal surface to a near-spherical shape, and 
the much slower timescales that involve a passage over energy 
barriers, such as the exchange of cell-cell neighbors. Such 
global topological rearrangements are associated with the lon- 
ger timescale T. 

Assessing the equilibrium or out-of-equilibrium nature of 
cystogenesis by means of a direct measurement of relaxation 
timescales is a daunting task. Here we propose a simpler method 
that allows unambiguous experimental discrimination between 
the two regimes. Assuming that the cyst is a monolayer and 
that spherical symmetry is preserved, the two opposite limits of 
equilibrium and strong deviation from equilibrium yield quite 
different predictions for the distribution of cell-cell contacts 
that characterizes the topology of the spherical monolayer. 

If mechanical equilibrium holds, the problem of energy 
minimization with given cell volumes can be mapped on the 
problem of tiling a sphere with polygons of a given surface. 



This problem has been solved by a combination of mathematical 
arguments and extensive computer searches (Cox and Flikkema, 
2010). For a cyst composed of several cells N > 13, the equilib- 
rium topology is given by the "soccer-ball" configuration made 
of 12 pentagons and N — 12 hexagons. 

In the opposite regimen mechanical equilibrium is never 
reached. In this case, for a 2D monolayer, our model predicts 
that the probability distribution of cell-cell contacts should re- 
duce to that found in Gibson et al. (2006) and Dubertret et al. 
(1998), as mechanical relaxations become negligible compared 
with cell division. Indeed, the models described in Gibson et al. 
(2006) and Dubertret et al. (1998) consider the topological dis- 
tribution of polygonal tilings emerging as a sole effect of random 
cell divisions, without taking into account mechanical relax- 
ation. This nonequilibrium topological distribution is consistent 
with the observed topologies for a large class of planar metazoan 
epithelia (Lewis, 1928; Rivier et al., 1995; Dubertret and Rivier, 
1997; Dubertret et al., 1998; Gibson et al, 2006; Li et al., 2012; 
Fig. SIN). The same arguments hold unchanged for curved epi- 
thelial structures, provided that the monolayer structure is pre- 
served and that the number of composing cells is large. 

Our numerical simulations qualitatively and quantitatively 
reproduce the theoretically predicted cyst geometries. Indeed, 
for x ^ 1, the nonequilibrium distribution emerges, as shown in 
Fig. 1 (A and B), and the corresponding topological distribution 
also corresponds to the predictions of Gibson et al. (2006) in this 
limit case (Fig. 1 D). In the opposite regimen (x 1), as shown 
in Fig. 1 C, the basal surface topology complies with the mini- 
mal energy tessellation of the sphere (Fig. 1 D, inset). 

Thus, the theoretical model predicts two opposite scenar- 
ios: if cysts grow under mechanical equilibrium, their surface 
topology has to be that of a soccer ball (Cox and Flikkema, 
2010). Conversely, if cysts grow out of mechanical equilibrium, 
their surface topology has to be the one predicted in Gibson 
et al. (2006). 

Cell-cell contact topology unveils 
the out-of-equilibrium nature of 
epithelial morphogenesis 

To test whether cysts grow under mechanical equilibrium or not, 
we segmented confocal sections of experimental MDCK cysts 
and reconstructed the complete 3D structure of several cysts (see 
Fig. 1, E and F). We measured the topology of the cell aggregate 
by computing the neighbor map on the outer surface of the cysts. 
As shown in Fig. 1 D, the topological distribution of MDCK 
cysts is very close to the one given in Gibson et al. (2006) and 
definitely different from the equilibrium distribution (Fig. 1 D, 
inset). Therefore, we conclude that in real cysts, cell division is 
rapid with respect to mechanical relaxation, i.e., real cells are 
growing in the x 1 regimen. Experimental cystogenesis oc- 
curs out of mechanical equilibrium. This theoretical conclusion 
is in agreement with existing data on real-time imaging of cyst 
growth showing that, both in cysts and in 2D epithelial layers, 
cell rearrangements (i.e., neighbor exchange) are rare (Gibson 
et al, 2006; Pearson and Hunter, 2007; Puliafito et al., 2012). 

The direct experimental measurement of relaxation times 
is practically challenging, as it is difficult to isolate relaxation 
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Figure 1 . Cysts are spherical monolayers out 
of mechanical equilibrium. (A) Computational 
simulation of cyst morphogenesis from 1 (top) 
to 32 cells (bottom), under nonequilibrium con- 
ditions and with control of the cell-division axis. 
The cross section (right) shows the presence of a 
single spherical lumen. The left column is the ex- 
ternal overview of the cyst. Equivalent time from 
seeding is around 4 d. (B) Numerical simulation 
of a cyst with 32 cells out of equilibrium. (C) The 
same cyst as in B, after relaxation to equilib- 
rium. (D) Analytical nonequilibrium topological 
distribution (orange; Gibson et al., 2006) and 
numerical simulations results (green). Normalized 
frequency of contacts for two MDCK cysts (n = 
53 cells, blue; and n = 1 1 1 cells, cyan) obtained 
by 3D reconstructions of confocal scans, (inset) 
Theoretical equilibrium topology distributions 
for two cysts of n = 53 (black) and n = 1 1 1 cells 
(gray). The result of equilibrium numerical simu- 
lations is identical (not depicted). (E) Four con- 
focal scans of the same cyst from bottom (left) 
to equator (right), stained with DAPI (blue) and 
phalloidin (green). Time from seeding: 4 d. Bar, 
10 urn. (F) 3D computer-assisted reconstruction 
of the cyst geometry of E from the full confocal 
scans. External overview (left) and cross sec- 
tion (right). Each cell is identified with a unique 
color for ease of visualization. Bar, 1 0 urn. 
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from other naturally occurring processes in the cyst, such as cell 
division or apoptosis. However, simulations allow us to easily 
perform such experiments numerically, and the details of out- 
of-equilibrium growth can be investigated. To this aim, we have 
simulated the growth of a cyst up to a size of a given number 
of cells, w C eiis, blocked cell division, and let the cyst relax me- 
chanically. Relaxation toward the equilibrium topology (Fig. 2, 
A and B) is a very slow process, as several energetically unfavor- 
able random cell movements have to take place one after another 
to overcome the many local energy barriers that separate the 
original nonequilibrium distribution from the globally energy- 
minimizing configuration (Fig. 2 A, insets). The time needed to 
reach the predicted equilibrium configuration increases exponen- 
tially with the number of cells in the cyst, as shown in Fig. 2 C. This 
finding confirms that cysts evolve in a rugged energy landscape 
and are geometrically frustrated systems (Mezard et al., 1987). 
Geometrical frustration means that local moves that would, by 
themselves, tend to minimize energy are actually forbidden or 
strongly discouraged by the constraints acting on the system. As 
a result, global rearrangements are required, with ensuing long 
times for relaxation to equilibrium. 

Quasi-equilibrium growth can be obtained if the cell di- 
vision time is significantly slowed down (Fig. 2 D), as this al- 
lows mechanical relaxation to take place between subsequent 
cell divisions. 

Although the typical nonequilibrium distribution described 
in Gibson et al. (2006) has been found in a variety of biological 
contexts, it has never been proved whether it is relevant in the 
context of developing epithelial structures in mammals. To test 
this hypothesis, we studied the topological distribution within 
tubular epithelia in developing mouse embryos. To this aim, we 
performed whole mount staining of lung and kidney samples 
and imaged them on a confocal microscope. After segmenting 
the data, the topology was calculated for rounded cyst-like 



alveoli from the lung (Fig. 3 A), the tubular ureter (Fig. 3 B), 
and the tips of the ureteric tree (Fig. 3 C) from the developing 
kidney. The analysis of the ureteric tree tips is a particularly rel- 
evant comparison to the in vitro cyst data, as MDCK cells were 
derived from an adult kidney. We found that the topological 
distributions from these tissues (Fig. 3 D) are indistinguishable 
from those found in an MDCK cyst (see Fig. 1 and Table S2). 
Therefore, the topological distribution is a universal fingerprint 
of nonequilibrium dynamics in cyst-like and tubular epithelial 
tissues both in vitro and in vivo. 

Relaxation of cortical tension in epithelial 
cysts favors equilibrium topology 

Our theory predicts that cysts could be brought to quasi-equilibrium 
conditions by slowing down proliferation, thus allowing more 
time for topological relaxation, or by favoring the passage over 
topological barriers that separate distinct configurations. 

The Rho-associated protein kinase (ROCK) inhibitor 
Y-27632 is known to relieve cortical tension (Matthews et al., 
2006) through inhibition of actomyosin contractility, and we 
therefore expect it to ease topological rearrangements and in- 
duce the appearance of equilibrium phenotypes. Indeed, theo- 
retically, ROCK inhibition would lower the overall mechanical 
relaxation time T, making x smaller. 

To verify the validity of our theoretical prediction, we 
compared the topology of growing cysts untreated and treated 
with ROCK inhibitor at different times from seeding, shown 
in Fig. 4 A. First, we evaluated the effect of ROCK inhibi- 
tion on cysts by using GFP-myosin light chain (GFP-MLC) 
as an indicator for cortical tension (Yu et al., 2003). As shown 
in Fig. 4 B, MLC-GFP is concentrated at the basal surface in 
untreated cysts (a quantification is given in Fig. S5). When 
Y-27632 is applied from the first day of the experiment, MLC 
is no longer accumulated at the basal surface, and instead is 
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Figure 2. The ratio of relaxation time to cell 
duplication time controls the approach to equi- 
librium. (A] Time trace of the energy during a 
numerical simulation of relaxation toward the 
equilibrium phenotype. The growth of a cyst 
with no control on spindle orientation is simu- 
lated up to 16 cells. Afterward, cell divisions 
are blocked and purely relaxational dynamics 
are generated. Note the presence of multiple 
sequential energy barriers during the conver- 
gence to the equilibrium phenotype. Time is 
measured in cell division time units. Note that 
although the single lumen is rapidly attained, 
the convergence to the full equilibrium topology 
takes much longer. (B) Topological distribution 
for the simulated relaxation of a 12-cell cyst. 
(C) The time T needed to reach the equilibrium 
configuration in the same in silico relaxation ex- 
periment as in A, starting from cysts composed 
of different numbers of cells, grows exponen- 
tially with the number of cells, (insets) Distribu- 
tion of cell-cell contacts at equilibrium for each 
value of n ce || s . (D) We simulated the growth of 
a normal cyst up to 16 cells (red trace and 
circles), and from that point on we let the cyst 
grow with different division rates. To compare 
different cysts with different division rates, we 
plot here the relative distance from equilibrium 
and nonequilibrium as a function of the number 
of cells. Slower cell divisions yield a topological 
shift toward equilibrium. 



diffuse in the cytoplasm, which is consistent with an alteration 
of cortical tension. 

To compare topologies in treated and untreated cysts, we 
digitally reconstructed >50 cysts, and, for each cyst, we mea- 
sured the topology of the outer surface and the number of com- 
posing cells. Representative images of experimental cysts are 
shown in Fig. 4 A and the complete set of reconstructed cysts is 
shown in Fig. S2. The comparison of theoretical predictions and 
experimental data for a set of cysts is summarized in Fig. 5 A. 
Cell-cell contact statistics of ROCK-inhibited cysts are sig- 
nificantly closer to the soccer-ball topology, which is the hall- 
mark of equilibrium, as shown in Fig. 5 (B-D). Conversely, as 



expected, the contact statistics of control cysts are closer to the 
out-of-equilibrium topology. 

To quantify the degree of equilibrium, we computed for 
each reconstructed cyst the normalized topological distribution 
and calculated the distance from equilibrium and nonequilib- 
rium distributions (see Materials and methods). The distances 
from equilibrium of cysts under different treatments at 1-5 d 
from seeding are shown in Fig. 5 C. 

According to our theoretical predictions, the time needed 
to reach mechanical equilibrium configuration depends on the 
number of cells in the cysts rather than on time (see Fig. 2 C). 
Thus, we expect the effect of ROCK inhibition to depend on the 
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Figure 3. Cystic and tubular epithelia in vivo and in vitro 
share the same topological distribution. Confocal optical sec- 
tions from fixed and stained whole-mount embryonic mouse 
kidney and lung tissues enable direct comparison of topolog- 
ical distributions. Representative mid, top, or bottom sections 
and the relative segmentation are shown for the three tissue 
types we analyzed. (A, left) Midsection of an El 3.5 lung 
alveolar tip stained with E-cadherin (green) and DAPI (blue). 
(A, right) Bottom section (top) and corresponding segmenta- 
tion (bottom) with E-cadherin staining only. Number of seg- 
mented cells in the image: 96. (B) Midsection of developing 
ureteric tubule in an El 3.5 mouse stained with E-cadherin 
(left), bottom section (middle), and corresponding segmenta- 
tion (right). Number of segmented cells in the image: 88. 
(C) Tip of a developing ureteric tree in the kidney of an 
El 3.5 mouse, stained with E-cadherin (green) and DAPI 
(blue; left), together with a top section (top right), and cor- 
responding segmentation (bottom right). Number of seg- 
mented cells in the image: 60. (D) Normalized topological 
distributions were reconstructed for three different lung tissue 
samples and four kidney tissue samples. All samples follow 
the nonequilibrium theoretical distribution. For the sake of 
comparison, the corresponding equilibrium topological dis- 
tribution is shown. Bars, 10 urn. 
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Figure 4. ROCK inhibition alters the cellular spatial dis- 
tribution of MLC. (A) Representative micrographs of top 
(top) and mid sections (bottom) of cysts fixed at days 3 
and 5 grown under no treatment (c3d and c5d) and 
treated with Y-27632 (Y3d and Y5d). Cysts were stained 
with DAPI and phalloidin. Bar, 10 urn. (B) We used GFP- 
MLC as an indicator for cortical tension (Yu et al., 2003). 
MDCK cells expressing GFP-MLC were cultured in Matri- 
gel solution. Aphidicolin and Y-27632 were applied sepa- 
rately from day 0 or day 3. Cells were then fixed at day 5 
and stained for PCX and nuclei. In control cysts, MLC is 
concentrated at the basal surface in an untreated cyst. 
When Y-27632 is applied from day 0, MLC is no longer 
accumulated at the basal surface, and instead is diffuse in 
the cytoplasm. In contrast, Aphidicolin does not change 
the localization pattern of MLC, which remains mainly at 
the basal surface. For a more detailed quantification of this 
effect see Fig. S5. Bar, 20 urn. 
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number of cells composing the cyst. Indeed, as shown in Fig. 5 
(B and C), Y-27632 is very effective on cysts composed by 
n < 40 cells, and its effect becomes weaker for larger aggre- 
gates. Note that for both treated and untreated samples, cysts 
containing a smaller number of cells are generally closer to 
equilibrium than bigger cysts, as expected. 

Treatment with Y-27632 has been reported to slow down 
cell divisions (Ishizaki et al., 2000). According to our theory, 
cell division slowdown can concur to relaxation toward the 
equilibrium configuration. To separate the contributions to this 
effect of cortical tension relief and cell division slowdown, we 
compared ROCK-inhibited cells with cells treated with the DNA 
synthesis inhibitor Aphidicolin. We performed experiments 
according to a variety of protocols, i.e., supplying the drugs in 
different combinations either from the first day of the experi- 
ment or from the third (Fig. 5 C). The observed slowing down 
of cell divisions induced by either Y-27632 or Aphidicolin is of 
a similar order (Fig. 5 E). However, topology is not significantly 
altered by Aphidicolin, as documented in Fig. 5 D. To explain 
this apparent discrepancy between our theoretical predictions 
and the experimental data, we recall the quantification of the 
effect of cell division slowdown on cell-contact topology (see 
Fig. 2 D). Indeed, the changes in the duplication time required 
to entail a significant effect on the relaxation toward equilibrium 
(on the order of 10-fold) are much larger than those experimen- 
tally accessible (for example, 30% to 40% in Fig. 5 E). Ac- 
cording to our theoretical model, the duplication time is indeed 



independent of the number of cells, whereas the relaxation time 
grows exponentially with the number of cells (Fig. 2 C). 

In summary, the observed equilibrium phenotype of 
Y-27632-treated cysts can most likely be attributed to the re- 
lease of cortical tension induced by ROCK inhibition, rather 
than to an effect on proliferation. 

Multi-luminal cysts result: from the 
combined effect of constrained 
out-of-equilibrium dynamics and 
disruption of mitotic spindle control 

The nonequilibrium character of cyst growth manifests itself 
with particular evidence in the dynamics of luminal spaces. It 
is well known that cysts grown with cells that are defective in 
establishing a proper direction for the cleavage plane are char- 
acterized by several, interspersed luminal spaces and the lack 
of a clear monolayer structure, as shown in Fig. 6 (A and B; 
Martfn-Belmonte et al., 2008). In vivo spherical epithelial 
structures equipped with such a phenotype would be severely 
impaired. A stringent control of the mitotic spindle orientation 
in the apico-basal direction therefore seems necessary. To test 
how sensitive the single-lumen phenotype is to perturbations 
in the alignment of the cleavage plane, we performed a large 
number of numerical simulations with decreasing control on 
the division geometry (Fig. 6, C-H). Indeed, our results indi- 
cate that larger errors in spindle positioning are associated with 
a dramatic increase of the emergence of aberrant phenotypes 
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Figure 5. Relaxation of cortical tension via inhibition of ROCK drives topology toward equilibrium, independent of cell division rate. (A) The table reports 
the experimentally observed and theoretically predicted topological configurations of a subset of control (top) and (bottom] MDCK cysts treated with 
Y-27632 to inhibit ROCK activity. Each line of the table represents, for a given number of cells in the cyst (top rows), the number of cells displaying four, 
five, six, seven, and eight neighbors. The topology of control cysts is never consistent with equilibrium and approaches the nonequilibrium topology for 
increasing numbers of cells, as expected from the theory. (B) After reconstructing the 3D structure of treated and untreated cysts, we calculated the distance 
from equilibrium and plotted it as a function of the number of cells. The topology of control cysts (closed blue triangles) is clearly closer to the nonequilibrium 
distribution, whereas Y-27632-treated cysts (closed red squares) are closer to equilibrium, especially when the number of cells is small. Here Y-27632 is 
given right after seeding. Open symbols are medians on series of data (n = 0-20, 20-40, 40-80, and 80-100), and lines are a guide to the eye. The 
data shown are from a total of n = 34 cysts, each containing between 10 and 90 cells, coming from a single experiment. (C) Topology of cysts at differ- 
ent times: untreated (top right), Y-27632 given from seeding (top left), Y-27632 given from day 3 (bottom left), and Aphidicolin given from day 3 (bottom 
right). Labels: c, control; Y, Y-27632; A, Aphidicolin; ld/2d/3d, 1/2/3 d; and c3dYld, untreated for 3 d, then treated with Y-27532 for 1 d; c3dA2d, 
untreated for 3 d, then treated with Aphidicolin for 2 d. Y-27632-treated cysts (top left) lie systematically in the half plane closer to equilibrium with respect 
to control cysts (top right) and contain on average a smaller number of cells (at a given time from seeding), due to the slowdown of cell division induced by 
Y-27632. Each point corresponds to the topology of a single cyst. (D) To account for the differences in topology, different groups of cysts were compared 
by means of a one-tailed Welch's test to assess for statistical significance, where the existence of a difference in the mean value of the topology was used 
as a null hypothesis. Cysts with 20-40 cells (left) and 40-60 cells (right) were grouped and compared. A statistically significant difference is seen between 
Y-27632-treated (left, n = 7; right, n = 6) and control (left, n = 7; right, n = 4) for any number of cells. As for other treatments, such as Aphidicolin (left, n = 5; 
right, n = 3) or Y-27632 given from day 3 (left, n = 1 ; right, n = 2), the topology is different from the Y-27632-treated cysts, resembling the untreated case, 
albeit not significantly. Asterisks indicate statistical significance: *, P < 0.05; **, P < 0.01 ; ***, P < 0.001 , etc. Error bars indicate standard deviations. 
(E) The effect of the different treatments on cell division was assessed by means of statistical comparison (Welch's Test) of different groups. The number of 
cells in the cysts at day 5 was compared in the different populations. All cases have n = 5 samples. Error bars indicate standard deviations. 
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Figure 6. The ratio of lumen coalescence to cell division times controls the monolumen or multilumen phenotypes. (A) We treated cysts with aPKC-PSi to 
induce polarity disruption and spatially disordered spindle positioning. Cells were seeded in growth medium supplemented with 50 pM aPKC-PSi, and fixed 
3 d later. Confocal scans for cysts treated with aPKC-PSi, which lack control over mitotic spindle positioning; stained for (3-catenin (green), labeling cell 
membranes; and ZO-1 (white) and PCX/gpl 35 (red), labeling luminal interfaces. Notice the multilumen morphology. Time from seeding: 3 d. Bar, 1 0 pm. 

(B) 3D reconstructions of the cyst in A from the full set of confocal sections. Several disconnected lumens are visible in the cross-section (right). Bar, 10 urn. 

(C) Mean number of luminal volumes and fraction of cysts with single lumen/multiple lumen phenotypes (inset), obtained from numerical simulations as a 
function of the tightness of the control over orientation of the cleavage plane, cp = 0°, mitotic planes always orthogonal to the luminal surface, cp = 90°, 
division planes chosen isotropically in three dimensions. Note that even for cp = 0°, it is still possible to get more than one lumen. Simulated cysts have 
32 cells. Error bars indicate standard deviations. (D) Numerical simulation of a multiple-lumen phenotype. The cleavage plane has been chosen at random, 
isotropically in 3D space. Equivalent time from seeding is around 5 d. (E) Mean fraction of cysts with single lumen/multiple lumen phenotypes (red/blue), 
obtained from numerical simulations as a function of the duplication time t, for different levels of the control over the orientation of the cleavage plane. 
The dynamics of several cysts were simulated for four rounds of division, and the fraction of normal cysts in the population was measured at the end 
of the simulation (16 cells), corresponding approximately to day 3 in the experiments. The duplication time is rescaled with to, which corresponds to 
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(Fig. 6 C, inset). To understand the dynamic origin of this phe- 
nomenon, we observed that for large proliferation rates, the 
process of lumen coalescence is bound to take place far from 
equilibrium, as new microlumens may be created at each cell 
division, and the number of new cell divisions per unit of time 
grows proportionally to the number of cell in the aggregate. In 
cysts with disrupted cell polarity, multiple lumens compete as 
aggregators of newly formed microlumens, as each new micro- 
lumen can merge with one of several different lumens in the 
cyst. If the typical time of lumen coalescence is small with re- 
spect to the cell division time, the monolumen phenotype is at- 
tained. Conversely, when coalescence is slow compared with 
cell division, the monolumen phenotype may never be reached. 
As normal MDCK cysts most often show a single lumen (Bryant 
et al., 2010), we expect the typical coalescence time to be faster 
than the cell division time. 

To further corroborate these observations we performed 
the following in silico experiment: we simulated the out-of- 
equilibrium growth of a cyst up to a given size, « C eiis, with dis- 
rupted orientation of mitotic spindle. The obtained multiluminal 
cyst was then relaxed, with cell divisions blocked. We mea- 
sured the time needed to reach the monolumen phenotype as 
shown in Fig. 6 F. Such times are faster than those required to 
reach global equilibrium (and in particular a soccer ball-like 
configuration), as shown in Fig. 2, which suggests the absence 
of significant energy barriers. The total time T m to reach the 
monolumen increases with the number of cells, as a larger num- 
ber of microlumens are generated. In our simulations, smaller 
lumens were observed to coalesce toward larger lumens by slid- 
ing along cell-cell interfaces, as depicted in Fig. 6 G. The slid- 
ing movement can take place without appreciably changing the 
global surface energy (Eq. 2), as the total area of cell-cell con- 
tacts does not change. We verified this prediction by simulating 
the coalescence of a single microlumen into a central lumen 
and by measuring the characteristic sliding time as a function of 
surface tension. Our results indicate that, indeed, coalescence 



time does not depend on surface tension (Fig. 6 H, inset). This 
is compatible with the experimental observation that ROCK 
inhibition does not have a significant effect on the speed of 
lumen coalescence (Fig. SI P). It is worth observing, however, 
that although the sliding of microlumens does not imply a sig- 
nificant energy expenditure, its final coalescence into a larger 
lumen provides an energy advantage, as a more symmetric and 
energetically favorable distribution of cellular interfaces is real- 
ized (Fig. 2 A, inset). The coalescence of lumens is therefore an 
energetically favored process. 

Theoretically, as in the case of the topology of cell con- 
tacts, slow cell division compared with lumen coalescence yields 
equilibrium growth, i.e., monoluminal phenotype, whereas fast 
cell division favors the accumulation of multiple lumens. To 
better characterize the role of the proliferation rate, we simu- 
lated the growth of a polarity-disrupted cyst with increasingly 
long division times, starting from a single cell. In normal cysts, 
lumen coalescence times are typically faster than cell division 
times, thus allowing single lumen configurations to be reached 
after each cell duplication. This is also observed in experimental 
cysts (Fig. 1, E and F). For cysts that have faster cell prolifera- 
tion, this ceases to be true, giving rise to multiple lumen pheno- 
types (Fig. 6 H). 

To experimentally validate our theoretical predictions, we 
quantitatively characterized the maturation and dynamics of lu- 
mens in MDCK cysts (Fig. SI P). In control cysts, the inci- 
dence of single lumens at day 2 is ^50%, while it reaches 80% by 
day 5. As monoluminal cysts become more frequent, the num- 
ber of immature and multiple lumens decreases. 

To verify whether lumen dynamics are actually independent 
of cortical tension, as suggested by our simulations (Fig. 6 H, 
inset), we generated multiluminal cysts by using an aPKC pseudo- 
substrate inhibitor (aPKC-PSi). This inhibitor contains a pseudo- 
substrate for atypical PKC that has been myristoylated to be cell 
permeable and alters the correct orientation of cell division (Qin 
et al., 2010). Cells were then treated with aPKC-PSi or Y-27632, 



normal MDCK cysts (20 degrees of error on spindle positioning and 30% of normal phenotype]. For t/t 0 > 1 (cell division slower than control), eventu- 
ally 100% of the cysts are normal, regardless of cell division orientation. For t/to, slightly smaller than unity, the emergence of aberrant phenotypes is 
strongly increased, even in the case of perfect spindle positioning. (F) We numerically simulated the relaxation of a multiluminal cyst grown with aberrant 
spindle positioning to the monolumen configuration, with blocked cell divisions. Here we plot the time needed to reach the monolumen rescaled by the 
division time as a function of the number of cells, n co u s . As expected, this increases with the size of the aggregates, as more microlumens are generated. 
(G) Cartoon illustrating how the process of lumen coalescence can take place with no alteration in the total extent of lateral surfaces, and therefore without 
changing the energy defined in Eqs. 1 and 2. (H) To show the effect of cell division on lumen coalescence in numerical simulations, we grew cysts starting 
from a single cell with different division rates. Cysts are grown with normal, i.e., correctly oriented, cell divisions. During time, cells divide and form micro- 
lumens that tend to coalesce. Faster cell division rates do not allow lumens to coalesce, independently of polarity. Error bars indicate standard deviations. 
(I) Aphidicolin treatment rescues the multilumen phenotype induced by aPKC perturbation. To check whether the predictions of our model on the slowdown 
of cell division were correct, we fixed cysts at day 4 and experimented with several treatments. Cysts were grown with aPKC-PSi from seeding (aPKC-PSi), 
with both aPKC-PSi and Aphidicolin (aPKC-PSi+Aph), or were treated with aPKC-PSi for 2 d, then washed and either left untreated (aPKC-PSi WO d2) or 
treated with Aphidicolin (aPKC-PSi WO+Aph]. Each percentage is obtained by means of the indicated number of independent experiments, for a total 
number of analyzed cysts reported in the top y axis. aPKC-PSi-treated cysts show a significant decrease in monoluminal cysts (aPKC-PSi [52.2 ±6.1] vs. 
control [73.4 ± 3.8]; p-value, 1 .5 x 10~ 12 ). Quadruplicate experiments of aPKC-PSi treatment and Aphidicolin show a significant shift toward the control 
situation, i.e., the number of monolumens increases (aPKC-PSi vs. aPKC-PSi+Aph [66.1 ± 1 .4]; p-value, 9.8 x 10~ 10 ). Washing out aPKC-PSi restores the 
normal level of mono- versus multiluminal cysts (aPKC-PSi vs. aPKC-PSi WO [67.0 ± 3.7]; p-value, 7.5 x 1 0~ 7 ]. Aphidicolin given for 2 d after the washout 
was similar to the washout alone (aPKC-PSi vs. aPKC-PSi WO+Aph [66.2 ± 4.7]; p-value, 1 .5 x 1 0~ 5 ). Statistical significance was assessed by means of 
a two-tailed Fisher test on 2 x 2 contingency tables containing mono- and multilumen counts and indicated treatments. Values are indicated as percentage 
means of the independent experiments ± standard deviation. N indicates total number of analyzed cysts and n indicates number of biological replicates; 
n = 2 + 2 indicates two biological replicates and two technical replicates. In the legend, "Mixed polarity" indicates cysts with gpl 35/PCX localized to the 
exterior surface of the cyst, where the cell membrane contacts the ECM. "Lumen" refers to cysts with centrally localized lumens, scored on the basis of 
the localization of gpl 35/PCX. "AMIS" refers to cysts that exhibit an accumulation of gpl 35/PCX in vesicular structures that underlie the membrane where 
the lumen will develop. 
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and the effect on lumen dynamics was quantified. As shown in 
Fig. SI P, the treatment with Y-27632 has no significant effect 
on the fraction of monolumen versus multilumen cysts. 

We then verified the effect of cell division on lumen dy- 
namics. To this aim, we measured the frequency of the monolu- 
minal phenotype in cysts with disrupted cell polarity where the 
proliferation rate was slowed down by Aphidicolin (Fig. 6 I). 
Interestingly, the slowing down of cell divisions induced a sig- 
nificant increase in the frequency of the monolumen phenotype 
(p-value: 10" 5 ). Although the control had 73.4 ± 3.8% monolu- 
mens, incubation with aPKC-PSi for 4 d decreased single lu- 
mens to 52.2 ± 6.1%. Slowing down of division by including 
Aphidicolin largely rescued the formation of single lumens, 
bringing it back to 66.1 ± 1.2%. 

In conclusion, slowing down cell division partially rescues 
the multilumen phenotype induced by misoriented cell division. 

Discussion 

Epithelial tissue growth is a complex process involving both 
biochemical signaling cascades and mechanical forces. The co- 
ordination between these two aspects is crucial for the correct 
development and maintenance of a well-ordered architecture, 
whereas aberrant phenotypes are the hallmark of life-threatening 
human diseases, such as adenocarcinomas and polycystic kidney 
disease. In this paper we focus on cystogenesis as an in vitro 
model of 3D epithelial growth and describe a novel theoretical 
framework that discriminates between the two radically differ- 
ent dynamical regimes observable when the cell aggregate grows 
either close to, or far away from, mechanical equilibrium. 

We show that cell-cell contact topology can efficiently 
discriminate between these two regimes. By comparing theoret- 
ical predictions for cell arrangements and the 3D reconstruction 
of experimental MDCK cysts, we demonstrate that cystogen- 
esis occurs out of mechanical equilibrium. 

The out-of-equilibrium condition can be rephrased in the 
observation that topological rearrangements in cysts involving the 
exchange of next-neighbors are rare events, occurring on time- 
scales much longer than the duplication time. This leads to a pecu- 
liar topological distribution of cell-to-cell contacts (Fig. 1 D) that 
is the true fingerprint of a strong departure from equilibrium. 

Neighbor rearrangements require that the cysts transiently 
assume energetically unfavorable configurations, for instance 
during the shrinking of a basolateral cell-cell contact area that 
is required in the process of neighbor exchange. Random cell 
motility favors the crossing of such topological energy barriers, 
therefore playing the role that thermal agitation has in statisti- 
cal physics systems. Indeed, induced activation of ERK1/2 in 
MCF10A acini results in increased cell motility and neighbor 
exchanges (Pearson and Hunter, 2007). 

Our experiments on cyst-like and tubular structures show 
clear hallmarks of nonequilibrium (Fig. 1) both in vivo and 
in vitro. These features are not shared by all epithelial tissues. For 
instance, Drosophila wing epithelium at pupal stage shows a nearly 
perfect hexagonal packing, corresponding to the equilibrium 
configuration for planar tilings (Classen et al., 2005). Similarly 
highly ordered, close-to-equilibrium structures are observed in the 



fovea, lens, and corneal tissues of mammals (Bhat, 2001; Jalbert 
et al., 2003; Hofer et al., 2005). The achievement of this geometri- 
cal order has been related to an increased ability in rearranging 
bonds between cells (Amonlirdviman et al., 2005; Classen et al, 
2005) induced by various mechanisms of active cell remodeling, 
such as increased trafficking of cell adhesion molecules (Classen 
et al., 2005) and anisotropic accumulation of myosin II (Rauzi 
et al., 2008). Cell division slowdown may also be at play in these 
tissues. Overall, these contributions act to favor barrier crossing 
and the approach to an equilibrium configuration. 

In MDCK cysts it has been reported that ROCK inhibition 
rescues the orientation of polarity in polarity-disrupted cells (Yu 
et al., 2008). Here we have shown that the inhibition of ROCK 
activity gives access to the equilibrium topological distribu- 
tion. ROCK inhibition is known to relieve cell cortical tension 
(Matthews et al., 2006), thus lowering the energy barriers that 
separate topologically different configurations and therefore 
favoring the approach to equilibrium. We have found that the 
multilumen phenotype induced by aPKC-PSi treatment cannot 
be rescued with a ROCK inhibitor, as shown in Fig. SIP, which 
suggests that lumen coalescence and that topological rearrange- 
ments of cell-cell contacts are unrelated processes. 

Another important feature that discriminates equilibrium 
versus nonequilibrium phenotypes is the presence of single or 
multiple lumens. Our results indicate that when proliferation rates 
are high with respect to the characteristic time for lumen coales- 
cence, multiple lumens may accumulate. Conversely, when pro- 
liferation rates are low, multiple lumens tend to coalesce into a 
single lumen. Additionally, when proper alignment of the mitotic 
spindle is disrupted by perturbing the normal establishment of cell 
polarity, aberrant multiple lumen cysts are even more frequent. 

In this case as well, nonequilibrium plays a determinant 
role: our results show that if cystogenesis were a pure equilib- 
rium process, mechanical relaxation would ultimately always 
lead to the single-lumen phenotype independently of cell polar- 
ity, as shown in Fig. 6 E. 

Our model explains why even in the case of normal cells, 
a significant fraction of cysts does not show a single lumen: 
indeed out-of-equilibrium growth can give rise to aberrant 
phenotypes even in presence of fairly accurate positioning of 
the spindle. 

Cell-cycle signaling pathways are often altered in cancers, 
and accelerated proliferation represents one of the most typi- 
cal features of carcinomas. In this context, the observation that 
control of cell duplication time alone can effectively have an 
impact on the properties of cyst growth is particularly relevant. 
Our simulations indicate that faster cell divisions would give 
rise to aberrant phenotypes even in the case where the geometry 
of cell division is perfectly controlled. These observations lead 
us to speculate that cancerous architectural patterns could also 
emerge due only to faster cell divisions, through a polarity- 
independent mechanism, thus calling for more detailed experi- 
mental investigations. 

The present findings raise the following question. Since 
equilibrium dynamics would, by themselves, ensure a proper 
cyst geometry without the need for tight control over the mitotic 
spindle orientation, why does cystogenesis not take place in the 
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equilibrium regimen? The answer probably lies in the fact that it 
is difficult to reconcile the contrasting requirements of fast growth 
and fast relaxation. Allowing sufficient time for mechanical re- 
laxation to equilibrium inevitably makes the overall endeavor of 
cyst construction much slower. To make matters worse, the num- 
ber of local energy barriers to be overcome before reaching equi- 
librium increases exponentially with the number of cells, with a 
consequent increase in the relaxation time. Eventually the cyst is 
trapped into an aberrant configuration. Control of the orientation 
of cell division circumvents this topological stalemate. 

This leads us to speculate that evolution has warranted 
fast and orderly cyst growth through an accurate control of cleav- 
age plane orientation, thereby remedying to the frustrating con- 
straints imposed by mechanics. 

Materials and methods 

Experimental procedures 

MDCK (subcloned from ATCC-CCL34 as described previously; Lipschutz 
et al., 2001) were cultured at 37°C, 5% CO2, in Eagle's minimal essential 
medium supplemented with 10% FBS, with the addition of L-Glutamine and 
antibiotics. To grow cysts, culture slides (Falcon 3541 08; BD) were precoated 
with a layer of basement membrane extracts (Matrigel; BD). Cells were 
seeded with a density of 50-1 00 cells/mm 2 and fed every other day in com- 
plete media supplemented with 2% Matrigel. aPKC activity was perturbed by 
the addition of 50 pM C-Zeta pseudosubstrate, myristoylated (aPKC-PSi; Bio- 
source) in growth media (Martin-Belmonte et al., 2007). ROCK activity was 
perturbed by the addition of 20 pM Y-27632 (EMD Millipore) to the me- 
dium. Aphidicolin was used at 1 pM. Cells were fixed and stained according 
to the protocol used in Yu et al. (2005). The primary antibodies used were: 
mouse anti-gpl 35/PCX (a gift from G. Ojakian, Department of Cell 
Biology, SUNY Downstate Medical Center, Brooklyn, NY), rabbit anti- 
B-catenin (Santa Cruz Biotechnology, Inc.), and rat anti-ZO-1 (R40.76; a 
gift from B. Stevenson, Nationwide Children's Hospital, Columbus, Ohio). 
Alexa fluorophore-conjugated secondary antibodies (Molecular Probes), 
Hoescht and DAPI (to counterstain nuclei), and Alexa Fluor-conjugated Phal- 
loidin were used. Whole-mount samples were dissected from embryonic day 
1 3.5 (El 3.5) C57/BI6 mice and fixed in 4% paraformaldehyde for 1 0 min. 
Samples were washed in phosphate-buffered saline with 0.1 % Triton X-100 
(PBTX), blocked (PBTX + 10% FCS), and then incubated for 48 h in primary, 
then secondary antibody solutions, washing with PBTX in between. Mouse 
anti-E-cadherin (No. 610181; BD), and rabbit anti-ZO-1 (No. 40-2300; 
Life Technologies) primary antibodies were used and detected with Alexa 
flurophore-conjugated secondary antibodies. Samples were then imaged on 
a laser scanning microscope (LSM510 and Leica SP2 equipped with PMT 
detectors, 40x oil objective lens, 1 .25 aperture; Carl Zeiss). The expression 
vector for MLC-GFP was transiently transfected into MDCK using Lipo- 
fectamine, according to the manufacturer's recommended procedure, and 
imaged after 24 h. 

Image analysis 

Image segmentation of cysts with formed lumen was performed with the help 
of Image], SciPy, Mayavi, and MATLAB software. To measure topological 
distributions in vivo, sections from confocal stacks were semiautomatically 
segmented, and neighbors were counted by means of MATLAB custom-writ- 
ten algorithms. Cysts were visualized by 3D rendering with a custom written 
code in SciPy. Luminal surfaces were measured from segmented stacks by fitting 
the data with triangulated surfaces and calculating their areas. In Figs. 2 D 
and 5 B, the distance from equilibrium was measured as cL/(d eq + d„ eq ), 
where d aq and d neq are the Euclidean distances of the vector x° (1 = 4, 5, ■ ■ •) 
of experimentally observed frequencies from the vectors x ne 9 and x ?9 of 
the theoretically predicted frequencies for nonequilibrium' (Gibson et al., 
2006) and equilibrium (Cox and Flikkema, 2010), respectively. 

Theoretical model 

We assume that the mechanical energy associated with a given cyst 
geometry is 

£ m = X ( aS cell-cell + P S cell-ECM + r S cell-l U men ) < (1) 
cells 



where S ce n_ CB n denotes the contact area between two cells (and similarly the 
contact area with ECM and lumen). The surface tensions a, fi, and y are 
hierarchically ordered: y > f} > a > 0. This ordering implies that cell-cell 
adhesion is energetically more favored than cell-ECM contact, leading 
to the formation of a spherical cell aggregate. The tendency to minimize 
mechanical surface energy is hindered by physical constraints, such as 
the limited variability of cell volume, and biochemical ones, such as the 
presence of signaling domains that control the size of the apical regions 
(Martin-Belmonte et al., 2007; Veglio et al., 2009; Bryant et al., 2010; 
see also Fig. SI). These two constraints are implemented through addi- 
tional terms C v (V ce ii), where V ce ii is the cell volume, and C A (S CIi |i_| umBn ). These 
penalization terms have a single minimum at a given target cell volume 
and target cell-lumen area, respectively. The total energy therefore reads 

E = E m + Z[C v (V cell ) + C A (S celWumen )]. (2) 

cells 

When a division occurs, a cell is split in two along a cleavage plane 
picked randomly in the plane parallel to the luminal surface. For cells lack- 
ing proper spindle positioning, the orientation of the cleavage plane is 
chosen completely (or at different extents) at random in 3D space. The 
model can be easily extended to the case of shape-dependent cell division 
(Gibson et al., 201 1), without substantially altering the predictions. After 
each division, the energy is redefined according to the new cell configura- 
tion. The constraint contribution to the energy has a much shorter relax- 
ation time than the surface energy part, so that during the dynamical 
evolution, cell volume and apical surfaces are only transiently different 
from their target values. 

Numerical simulations 

Cystogenesis was simulated using a dynamic Potts model (Potts, 1952; 
Wu, 1982; Graner and Glazier, 1992; Glazier and Graner, 1993). On 
a 3D cubic lattice, each site is occupied by a "color" variable that takes on 
integer values. To each value there corresponds an individual cell. At each 
division a new color is added. The exchange energies on neighboring sites 
provide a lattice approximation of the surface energies of cell-cell, cell- 
lumen, and cell-ECM contacts described in Eq. 1 . Target cell volume and 
luminal surface are introduced on the basis of experimental observation 
and enforced by additional terms in the energy. The target cell volume and 
lumen areas are set to V 0 = 1 ,200 urn 3 and S 0 = 80 urn 2 , according to ex- 
perimental values (Fig. SI). Simulations always started from a single cell. 
Duplication times are chosen from a uniform distribution with mean t and 
coefficient of variation 0.1. To simulate normal cells, the division axis is 
chosen at random in the plane tangent to the lumen surface, whereas in the 
aberrant case it is allowed to depart from the tangent plane to various de- 
grees. We verified that our results do not depend significantly on the pre- 
cise values of the surface tension parameters in Eq. 1 , as long as the 
hierarchy y > (J > a is respected. The results are similarly independent of 
the precise form of the constraints in Eq. 2, as long as they have single 
minima and a relaxation time much shorter than the typical relaxation time 
of the surface energy term. 

Stochastic algorithm 

The dynamic Potts model is driven by a discretization of the energy as- 
signed by Eqs. 1 and 2. We consider a 3D cubic lattice of N = n 3 sites, 
with n = 100-300 and periodic boundary conditions. We assign to each 
site f a variable S/ whose value indicates if the site is occupied by ECM (s, = 0), 
lumen (s, = — 1), or one of a collection of n c cells (/=!,..., n c ). Every 
color cluster corresponds to a cell. We refer to the s, values as "spins" or 
"colors," following the traditional nomenclature of statistical physics. Dis- 
cretized surface energies are computed by the formula 

N r 
H m = Z4,ct ) . 1_S cr i ,<T j < 

where the sum runs over the edges that connect each site / to its 26 cubic 
neighbors. The symmetric tensor -'0,0, describes the energy cost of having 
two neighboring sites occupied by different colors. It takes on different 
values depending on whether an edge joins two cells, a cell and the ECM, 
or a cell and the lumen. There is no cost for edges joining two sites with 
the same color. To reduce the anisotropy induced by the cubic lattice, 
the interactions are weighted depending on the edge length. The closest, 
next-closest, and next-to-next-closest neighbors have weights 1 , 1 /I, and 
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0, respectively. Target lumen area is enforced by adding to H m the penal- 
ization term 



^aZ(A-s 0 )^, 

1=1 

where the Lagrange multiplier A A satisfies a single requirement: it must be 
large enough to provide a relaxation time for the constraint contribution 
much faster than the characteristic time for the relaxation of the surface en- 
ergy. Target volumes are enforced by the Monte Carlo dynamics by allowing 
only moves that tend to decrease the distance between the actual cell volume 
and the target volume Vb, as described in Creutz (1 983]. We start the simu- 
lations from a single spherical cell with "color" 1 . Spins relax according to 
a Monte Carlo dynamics with inverse temperature (3 (Landau and Binder, 
2000). The precise value of (3 is not relevant as long as (3 is large enough to 
smear lattice anisotropies associated with the cubic discretization and small 
enough to allow the formation of stable spin clusters. At each time step, 
a randomly chosen spin is "flipped" to another color chosen among those 
taken by its neighbors, corresponding to an elementary extension or retrac- 
tion of one cell to the expenses of its neighbors. This random move is ac- 
cepted if the change in energy AH is negative. It is rejected with probability 
1 — e~ p4H if AH > 0. When a move is accepted, another spin is flipped, and 
the procedure is repeated iteratively until the division time, chosen for each 
cell from a uniform probability distribution with mean t and coefficient of 
variation 0.1, is reached. Then, the original cell is allowed to grow to the 
double of the target volume V 0 and split in two along a randomly chosen 
cleavage plane. Two new colors are assigned to the newly formed cells. A 
germ of luminal surface is positioned at the interface between the two cells 
by assigning —1 spins in a localized circle around the center of the cleav- 
age plane. The orientation of the cleavage plane is always randomly cho- 
sen, but the plane is constrained to be perpendicular to the cell apical 
surface for wild-type cells, and is instead completely unconstrained for cells 
lacking any control on mitotic spindle orientation. To determine the amount 
of error that is necessary to destroy the monolumen architecture, the interme- 
diate case was also considered, where the positioning of the cleavage 
plane perpendicular to the apical surface could be affected by an angular 
error of increasing amplitude. 

Solvable model for cyst radii 

The configuration of a cyst of N cells, subject to the constraints that cells 
have volume Vand apical surface a can be calculated analytically in the 
limiting case of sharp constraints. Indeed, the internal and external radii 
(rand R, respectively) obey: (4it/3)(P 3 - r 3 ) = NV, Am 2 = Na. From these 
two relations we obtain r = JNa / 4ti and R = [(N o/4ir) 3/2 + 
(3/4i7)NVV /3 . 

Image analysis 

DAPI- and Phalloidin-stained cysts stacks were acquired by confocal mi- 
croscopy, with LSM510 (Carl Zeiss) and SP2 (Leica) microscopes (Fig. S3, 
A and D). The boundaries of cells in each slice were drawn semiauto- 
matically by looking at the original images. Both nuclei position and actin 
staining were used to this aim. Each cell was tracked across the stack by 
assigning it a unique color in the stack. All cells in the cysts were then 
uniquely identified, as shown in Fig. S3 (B and E). Black traits correspond- 
ing to cell boundaries were then fused to neighboring cells by a standard 
erosion algorithm. A median filter was applied to remove isolated pixels 
and remaining segmentation errors (Figs. S3, C and F). Slices were then 
mounted together in a 3D image at the correct spatial resolution (i.e., ac- 
counting for different resolutions in the x-y plane and in the z direction). 
The result was then visualized by a rendering routine written in Python 
and MayaVi (Figs. S3, G and H). To extract the outer (or inner) topology, 
a 3-pixel-thick shell was extracted and contacts were then computed on 
the shell by calculating the neighbor map. Results were then carefully cor- 
rected by eye on the full 3D structure to identify possible errors coming from 
the finite spatial resolution. To quantitatively characterize MLC-GFP spatial 
distribution, we segmented midsections of digital images of stained cysts 
by using DAPI, Phalloidin, and GP135 staining. The three channels were 
used to reconstruct a binary mask identifying the section of the cyst. Within 
the mask, we identified cells by segmenting nuclei and then applying the 
watershed transform into the mask. Two circular layers spanning a quarter 
of the thickness of the cyst (one starting from the basal and the other starting 
from the apical side) were digitally constructed. The ratio of fluorescence 
intensity of the MLC-GFP channel between the basal and apical quarter was 
then used as a quantification (Fig. S5). 



Statistical analysis of topological distributions 

The statistical likelihood that a given experimental distribution originates 
from nonequilibrium dynamics is determined as follows. The probability of 
observing n, counts in a histogram for polygons of / sides, with theoretically 
expected frequencies \l„ is multinomial: 



P(n ] ,...,n k lw,...,Hi) = jp^I>/' 



(3) 



and is zero whenever one of the jj./ is zero. Thus, we have to restrict the 
comparison to nonzero entries of the expected distribution, i.e., to poly- 
gons with 4, . . . , 9 sides. We refer here to one of the distributions re- 
ported in Gibson et al. (201 1 ; Fig. 6 B, gray and blue lines). In particular, 
we do not take into account refinements of that model (such as, for exam- 
ple, those described in Gibson et al. [2006] or Sandersius et al. [201 1]). 
The original Gibson distribution has been derived by considering random 
unbiased cell division, which is known to be an oversimplification. How- 
ever, the effect of all the corrections that have been made afterward is to 
change the final distribution by a few percentage points, and does not rep- 
resent, from a quantitative point of view, a major difference. In our work 
we do not discuss any of these factors, such as size-dependent cell division 
rates, or shape-dependent control of mitotic spindle. The two opposite 
cases treated in the present work are the equilibrium and out-of-equilibrium 
cases, which do show significantly different distributions. The likelihood of 
observing a given realization, given the theoretical distribution, and the 
corresponding errors can be evaluated by Bayes' formula: P(|xln) = 
P(nl|x)P 0 ((j.)/P(n). Following a maximum likelihood approach, we assume a 
uniform prior, i.e., Po = const, and maximize P(|xln) by maximizing P(nljju). 
Corresponding values of |jl * such that the likelihood function is maximized 
are those satisfying: 



logP + X 1-Xn, 



and 



\ogP + X l-2>f 



= 0 



= 0, 



where A is a Lagrange multiplier used to impose Z|x; = 1 . The solution is 
|x* = n,/N. We estimate the Euclidean distance between the estimated |jl * 
and the expected p. as 



To estimate the corresponding error one has to calculate the Hessian: 
<) 2 logP/()n,,fl(j,j = (N 2 /n,)S, ( . Here, the curvature of the Hessian represents 
the error. For the k — th bin, this is -Jr\ / N . Thus the error on the dis- 
tance is 



Ad = d 



For the observation to be consistent with nonequilibrium dynamics, d must 
be compatible with zero with a standardized test, i.e., d ± 3Ac/ must in- 
clude 0 at the 99% confidence level. To implement these quantifications, 
we generated several synthetic distributions by extracting random numbers 
according to the probabilities given in Table SI. The corresponding dis- 
tances and probability values are given in Fig. S4, C and D. These statisti- 
cal methods, applied to the data extracted from Figs. 1 and 3, give the 
results reported in Table S2. The same methods, applied to the data ex- 
tracted from Fig. 4 C, give the results reported in Fig. S4 (A and B). 

Online supplemental material 

Fig. SI provides quantitative data about measured epithelial geometry and 
topology. Fig. S2 shows digital 3D reconstructions of fixed and stained 
MDCK cysts. Fig. S3 details the image segmentation process. Fig. S4 con- 
tains data about the statistical analysis of cyst topology. Fig. S5 shows a 
quantification of GFP-MLC distribution in untreated and ROCK-inhibited 
cysts. Table SI contains predicted nonequilibrium probabilities for the 
number of cell-cell contacts. Table S2 shows distances of experimental 
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data from the nonequilibrium distribution. The following custom written rou- 
tines for 3D stochastic simulations and image visualization are provided in 
a zip file: potts3d.f is a Fortran code for dynamical simulations of 3D Potts 
model; image_analysis.py is a Python code for image analysis using SciPy 
and MayaVi libraries; invivo_topology.py is a MATLAB code for the analy- 
sis of in vivo data; and palette is the list of colors used by the routines in 
"image_analysis.py." Online supplemental material is available at http:// 
www.icb.org/cgi/content/full/jcb.201 305044/DC1 . 
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